Unsteady non-Newtonian hydrodynamics in granular gases 
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The temporal evolution of a dilute granular gas, both in a compressible flow (uniform longitudi- 
nal flow) and in an incompressible flow (uniform shear flow), is investigated by means of the direct 
simulation Monte Carlo method to solve the Boltzmann equation. Emphasis is laid on the identifi- 
cation of a first "kinetic" stage (where the physical properties are strongly dependent on the initial 
state) subsequently followed by an unsteady "hydrodynamic" stage (where the momentum fluxes 
are well-defined non-Newtonian functions of the rate of strain). The simulation data are seen to 
support this two-stage scenario. Furthermore, the rheological functions obtained from simulation 
are well described by an approximate analytical solution of a model kinetic equation. 
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I. INTRODUCTION AND STATEMENT OF 
THE PROBLEM 

Granular gases are dilute systems made of inelastic 
particles that can be maintained in a fluidized state by 
the application of external drivings to compensate for the 
dissipation of kinetic energy due to collisions. These sys- 
tems are always out of equilibrium and exhibit a wealth 
of intriguing complex phenomena [1—12]. They are im- 
portant from an applied point of view but also at the 
level of fundamental physics. As Kadanoff stated in his 
review paper [4], "one might even say that the study of 
granular materials gives one a chance to reinvent statis- 
tical mechanics in a new context." 

One of the most controversial issues in granular flu- 
ids refers to the validity of a hydrodynamic description 
[4, 13]. In conventional fluids, the densities of the con- 
served quantities (mass, momentum, and energy) satisfy 
formally exact balance (or continuity) equations involv- 
ing the divergence of the associated fluxes. In the case 
of granular fluids, however, energy is dissipated on col- 
lisions and this gives rise to a sink term in the energy 
balance equation. As a consequence, except perhaps in 
quasielastic situations, the role of the energy density (or, 
equivalently, of the granular temperature) as a hydro- 
dynamic variable is not evident. Both for conventional 
and granular fluids, the mass, momentum, and energy 
balance equations do not form a closed set due to the ap- 
pearance of the momentum and energy fluxes (plus the 
energy sink in the granular case) . On the other hand, by 
assuming "hydrodynamic" conditions, the balance equa- 
tions are closed by the addition of approximate constitu- 
tive equations relating the momentum and energy fluxes 



/[r,v,t|/ ] /[v|n,u,T] 



/o(r,v) 



+ 



time 



Kinetic Hydrodynamic 
stage stage 



* aavivas@unex.es; http://www.unex.es/eweb/fisteor/antonio_ 
astillero/ 

i andrcs@unex.es; http://www.unex.es/eweb/fisteor/andres 



FIG. 1. (Color online) Schematic description of the two-stage 
evolution of the velocity distribution function in a conven- 
tional gas. 



(again plus the energy sink in the granular case) to the 
mass, momentum, and energy fields. 

The simplest constitutive equations consist of replacing 
the fluxes by their local equilibrium forms, thus neglect- 
ing the influence of the hydrodynamic gradients. This 
gives rise to the Euler hydrodynamic equations, which 
fail to account for irreversible effects, even in the case 
of conventional fluids. This is corrected by the Navier- 
Stokes (NS) constitutive equations, where the fluxes are 
assumed to be linear in the hydrodynamic gradients. On 
the other hand, if the gradients are not weak enough 
(i.e., if the Knudsen number is not small enough), the 
NS equations are insufficient and, thus, nonlinear (i.e., 
non-Newtonian) constitutive equations are needed in a 
hydrodynamic description [14, 15]. 

In conventional fluids, applicability of a hydrodynamic 
description (Euler, NS, or non-Newtonian) requires two 
basic conditions, one spatial and another one temporal. 
On the one hand, one must focus on the bulk region of the 
system, i.e., outside the boundary layers, whose width is 
on the order of the mean free path. On the other hand, 
one must let the system age beyond the initial layer, 
whose duration is on the order of the mean free time. 
Let us consider the latter condition in more detail. In a 
conventional gas, the typical evolution scenario starting 
from an arbitrary initial state represented by an arbitrary 
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FIG. 2. (Color online) Sketch of the USF. 
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FIG. 3. (Color online) Sketch of the ULF for (a) a(t) > and 
(b) a(t) < 0. 



initial velocity distribution /o(r,v) proceeds along two 
successive stages [16]. First, during the so-called kinetic 
stage, the velocity distribution /(r, v,i), which depends 
functionally on /o, experiences a fast relaxation (lasting a 
few collision times) toward a "normal" form /[v|n, u, T], 
where all the spatial and temporal dependence occurs 
through a functional dependence on the hydrodynamic 
fields (number density n, flow velocity u, and tempera- 
ture T). Next, during the hydrodynamic stage, a slower 
evolution of the hydrodynamic fields takes place until ei- 
ther equilibrium or an externally imposed nonequilibrium 
steady state is eventually reached. While the first stage is 
very sensitive to the initial preparation of the system, the 
details of the initial state are practically "forgotten" in 
the hydrodynamic regime. Figure 1 depicts a schematic 
summary of this two-stage evolution in a conventional 
gas. 

The absence of energy conservation in granular fluids 
sheds some reasonable doubts on the applicability of the 
above scenario beyond the quasielastic regime. While the 
usefulness of a non- Newtonian hydrodynamic description 
in steady states has been validated by computer simu- 
lations [17-20], it is not obvious that a hydrodynamic 
treatment holds as well during the transient regime to- 
ward the steady state. 

In order to address the problem described in the pre- 
ceding paragraph, it seems convenient to focus on cer- 
tain prototypical classes of flows. Let us assume a (d- 
dimensional) granular gas with uniform density n(t), uni- 



form temperature T(t), and a flow velocity along a given 
axis (say x) with a linear spatial variation with respect 
to a certain Cartesian coordinate £, i.e., 



VjW;(r,i) 



a(t)5 ix 5j£. 



(1.1) 



Here a(t) is a uniform rate of strain. Two distinct possi- 
bilities arise: either I ^ x (say I = y) or £ = x. The 
first case defines an incompressible flow (V • u = 0) 
commonly known as simple or uniform shear flow (USF) 
[1, 11, 15, 17, 21-60], the associated rate of strain a being 
the shear rate. The second case is an example of com- 
pressible flow (V • u = a 7^ 0) that will be referred to as 
uniform longitudinal flow (ULF) [60-67] ; the correspond- 
ing rate of strain in this case will be called longitudinal 
rate. These two states are particular cases of a more 
general class of homo-energetic affinc flows characterized 
by d 2 Ui/dxjdx k = [21]. The USF and ULF flows are 
sketched in Figs. 2 and 3, respectively. 

Assuming that the velocity distribution function 
f(r,v,t) depends on the spatial variable £ only, the Boltz- 
mann equation reads 



Tt +Vt m =J[Lfl 



(1.2) 



where J[f, f] is the (inelastic) Boltzmann collision op- 
erator, whose explicit form can be found, for instance, 
in Rcfs. [68-70]. Multiplying both sides of Eq. (1.2) by 
{l,v, u 2 }, and integrating over velocity, we get the bal- 
ance equations for mass, momentum, and energy densi- 
ties, 



D t n = 



D t Ui 



dn 



1 dP e , 



di ' 



(1.3) 



(1.4) 



dui uq £ \ 
P ^ + ^'- (L5) 



In Eqs. (1.3)-(1.5), D t = dt+uede is the material deriva- 
tive, and the number density n, flow velocity u, tem- 
perature T, pressure tensor P^ , heat flux vector q, and 
cooling rate £ are defined by 



n(£, t) 



dv/(4v,f), 



(1.6) 



n(£,t)u(£,t) = / dvv/(4v,t), (1.7) 



n(£, t)T{l, t) = p{£, t) = ^trP(£, t), (1. 



P tJ (£,t)=m / &v[vi - Ui^t)}^ - Uj (l,t)]f{l,v,t), 



(1.9) 



3 



q(l,t) = ^ J dv [v-u(£,t)] 2 [v-u(l,t)}f(£, V ,t), 

(1.10) 



n(£,t)T(£,t)((e,t) 
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dv« 2 J[/,/]. (1.11) 



The first equality of Eq. (1.8) defines the hydrostatic 
pressure p, which is just given by the ideal-gas law in 
the Boltzmann limit. 

As said above, the density n(t) and temperature T(t), 
and so the hydrostatic pressure p(t), are uniform in the 
(fully developed) USF and ULF. On physical grounds, it 
can also be assumed that the whole pressure tensor (t) 
is uniform as well. Moreover, in the absence of thermal 
and density gradients, the heat flux can be expected to 
vanish. Taking all of this into account, as well as Eq. 
(1.1), the balance equations [Eqs. (1.3)— (1.5)] become 

n{t) = -n(t)a(t)5 xe , (1.12) 
a(t) = -a 2 (t)S xe , (1.13) 



f(t) = -^P x t(t)-t(t)T(t). (1.14) 

In the case of the USF (t = y), Eqs. (1.12) and (1.13) 
imply that both the density and the shear rate are con- 
stant quantities. As for the temperature, it evolves in 
time subject to two competing effects: viscous heating 
[represented by the term —(2/dn)aP xy > 0] and inelas- 
tic cooling [represented by £T]. Both effects eventually 
cancel each other in the steady state. 

In the case of the ULF (£ = x), the solution to Eqs. 
(1.12) and (1.13) is 



"0 



n(t) n Q 

~ TT = — , a(t) = 
a(t) ciq 1 + dot 



(1.15) 



where no is the initial density and ao is the initial longi- 
tudinal rate. In contrast to the USF, the sign of a{t) (or, 
equivalent ly, the sign of ao) plays a relevant role and de- 
fines two separate situations (sec Fig. 3). The case a > 
corresponds to a progressively slower expansion of the 
gas from the plane x = into all of space. On the other 
hand, the case a < corresponds to a progressively faster 
compression of the gas toward the plane x = 0. The 
latter takes place over a finite time period t = |ao| _1 . 
However, since the collision frequency rapidly increases 
with time, the finite period t = |ao| _1 comprises an in- 
finite number of collisions per particle [67]. Given that 
Pxx > 0, the energy balance equation [see Eq. (1.14)] im- 
plies that the temperature monotonically decreases with 
time in the ULF with a > 0. On the other hand, if a < 0, 
we have again a competition between viscous heating and 
inelastic cooling, so the temperature either increases or 
decreases (depending on the initial state) until a steady 
state is eventually reached. The main characteristic fea- 
tures of the USF and the ULF are summarized in Table 
I. 



TABLE I. Main characteristic features of the USF and the 
ULF. 

ULF (a > 0) ULF (q < 0) 



USF 


Inelastic cooling 


Yes 


Viscous heating 


Yes 


T(t)4.&|a*(t)|t 


Yes, if 




f ^ S\aP xy 
^ dp 


T(i) t & \a*(t)\l 


Yes, if 




^ dp 


Steady state 


Yes 



Yes 
No 
Yes 

No 

No 



Yes 
Yes 
Yes, if 



2\a\P x 
dp 

Yes, if 

2\a\P x 



dp 



Yes 



Regardless of whether the rate of strain a is constant 
(USF) or changes with time (ULF), the relevant param- 
eter is the ratio 



a (t) = -77T 



(1.16) 



between a(t) and a characteristic collision frequency 
v(t) oc n(t)[T(i)]V2, Note that the absolute value of a* (t) 
represents the Knudsen number of the problem, i.e., the 
ratio between the mean free path and the characteris- 
tic length associated with the velocity gradient [59, 60]. 
Since a(i)/n(t) = const both in the USF and the ULF, 
we have a*(t) oc [T(t)] _1 / 2 . Consequently, the qualita- 
tive behavior of |a*(t)| is the opposite to that of T(t), as 
indicated in Table I. 

The scenario depicted in Fig. 1, if applicable to a gran- 
ular gas in the USF or in the ULF, means that, af- 
ter the kinetic stage, the velocity distribution function 
/[r, v, t|/o] should adopt a hydrodynamic (or normal) 
form 



/[r,v,t|/ ] 
where 



n(t) 



2T(t) 



d/2 



f*(C(r,t),a*(t)), 

(1.17) 



C(r,t) 



v - u(r, t) 

V 2T ( t )/ m 



(1.18) 



is the peculiar velocity scaled with the thermal speed. For 
a given value of the coefficient of restitution, the scaled 
velocity distribution function /*(C, a*) must be indepen- 
dent of the details of the initial state fo and depend on 
the applied shear or longitudinal rate a(t) through the 
reduced scaled quantity a* only. In other words, if a hy- 
drodynamic description is possible, the form (1.17) must 
"attract" the manifold of solutions /[r,v,£|/o] to the 
Boltzmann equation (1.2) for sufficiently long times, even 
before the steady state (if it exists) is reached. Equation 
(1.17) has its counterpart at the level of the velocity mo- 
ments. In particular, the pressure tensor fy[i|/o] would 
become 

P t] [t\h]^n(t)T{t)P: 3 {a*{t)) (1.19) 
with well-defined hydrodynamic functions P*j(a*). 
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A few years ago we reported a preliminary study [58] 
where the validity of the unsteady hydrodynamic forms 
(1.17) and (1.19) for the USF was confirmed by means 
of the direct simulation Monte Carlo (DSMC) method to 
solve the Boltzmann equation and by a simple rheological 
model. The aim of the present paper is twofold. On the 
one hand, we want to revisit the USF case by present- 
ing a more extensive and efficient set of simulations, by 
providing a detailed derivation of the rheological model 
(which was just written down without derivation in Ref. 
[58]), and by including in the analysis the second visco- 
metric function (which was omitted in Ref. [58]). On the 
other hand, we perform a similar analysis (both compu- 
tational and theoretical) in the case of the ULF. This 
second study is relevant because, despite the apparent 
similarity between the USF and the ULF, the latter dif- 
fers from the former in that it is compressible, the two 
signs of a physically differ, and a steady state is possible 
only for negative a. 

The remainder of the paper is organized as follows. 
The formal kinetic theory description for both types of 
flow is presented in Sec. II within a unified framework. 
Next, Sec. Ill offers a more specific treatment based on a 
simple model kinetic equation. In particular, a fully ana- 
lytical rheological model is derived. Section IV describes 
the simulation method employed to solve the Boltzmann 
equation and the classes of initial conditions considered. 
The most relevant part of the paper is contained in Sec. 
V, where the results obtained from the simulations are 
presented and discussed. It is found that the scenario 
depicted by Fig. 1 and Eqs. (1.17) and (1.19) is strongly 
supported by the simulations. Moreover, the simple ana- 
lytical rheological model is seen to agree quite well with 
the simulation results. The paper is closed in Sec. VI 
with a summary and conclusions. 



A. Changes of variables 

We start by defining scaled time and spatial variables 
t and £ as 



where 



6{t) 



t = 6(t), e=9(t)l, (2.1) 
USF (£ = y), 



6{t) 



- l \n{l + aot), ULF (£ = £), 



USF (£ = y), 
ULF (I = x). 




Also, a new velocity variable v is defined as 
v = v - a(t)£e x , 



(2.2) 



(2.3) 



(2.4) 



where is the unit vector in the x direction and we 
have taken into account that a(t) = ao9(t) both in the 
USF and in the ULF. The velocity distribution function 
corresponding to the variables £, v, and t is 



f(£, v,t) = J-/(f )V) t). 
0{t) 



Consequently, 
1 df df 



- a S x 



\ d£ dv x 



9 2 dt dt 
where we have taken into account that 
6(t) = -a [6{t)] 2 5 xi . 

Similarly, 



II. BOLTZMANN EQUATION FOR USF AND 
ULF 



l_5/_9/_ df 



(2.5) 



(2.6) 



(2.7) 



(2.8) 



Let us consider a granular gas modeled as a system of 
smooth inelastic hard spheres (of mass m, diameter <r, 
and constant coefficient of normal restitution a), subject 
to the USF or to the ULF sketched in Figs. 2 and 3, re- 
spectively. In the dilute regime, the velocity distribution 
function f(£,v,t) obeys the Boltzmann equation (1.2). 
As is well known, the adequate boundary conditions for 
the USF are Lees-Edwards's boundary conditions [71], 
which are not but periodic boundary conditions in the co- 
moving Lagrangian frame [72] . The appropriate bound- 
ary conditions for the ULF are much less obvious. In 
order to construct them, it is convenient to perform a 
series of mathematical changes of variables. Also, we 
will proceed by encompassing the ULF and the USF in 
a common framework. 



j 2 J[fJ] = J[f,f]. (2.9) 

Inserting Eqs. (2.6), (2.8), and (2.9) into Eq. (1.2), and 
taking into account Eq. (2.4), one finally gets 

l + ^|- a °^(^> J[/ ~ /] - (2 ' 10) 

It is important to remark that no assumption has been 
made. Therefore, Eqs. (1.2) and (2.10) are mathemat- 
ically equivalent, so any solution to Eq. (1.2) can be 
mapped onto a solution to Eq. (2.10) and vice versa. 
While Eq. (1.2) describes a gas in the absence of external 
forces, Eq. (2.10) describes a gas under the influence of 
a nonconservative external force F = —maoV{e x . Note 
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that in the case of the ULF with ao < the finite time 
interval < t < |ao| _1 translates into the infinite scaled 
time interval < t^< oo. 

The density n(£,t), flow velocity u(£,t), temperature 
T(£,t), and pressure tensor Pij(£,t) associated with the 
scaled distribution (2.5) are defined analogously to Eqs. 
(1.6)-(1.9). The quantities with and without tilde are 
related by 



n(£,t) 



n{£,t) 



Pg(t,t) 
8(t) 



(2.11) 



u(£,t) = u(£,t) -a(t)£e x , T(£,t) =T(£,t). (2.12) 

At a microscopic level, we now define the USF (£ = y) 
and the ULF (£ = x) as spatially uniform solutions to 
Eq. (2.10), i.e., 



f(£,v,t) = f(v,t). 



(2.13) 



Thus, conservation of mass and momentum implies that 
n = const and u = const. Without loss of generality 
we can take u = 0. It seems quite natural that periodic 
boundary conditions (at £ = ±L/2) are the appropri- 
ate ones to complement the (scaled) Boltzmann equation 
(2.10) in order to ensure the consistency with uniform so- 
lutions (2.13), i.e., 



f(-L/2,v,t) = f(L/2,v,t). 



(2.14) 



Assuming uniform solutions of Eq. (2.10) and going 
back to the original variables, Eqs. (2.11) and (2.12) 
yield n(t) = 6(t)n, u(£,t) = a(t)£e x , T(t) = f(t), and 
Pi jit) = 6{t)Pij(t). Therefore, uniform solutions to Eq. 
(2.10) map onto USF (£ = y) or ULF {£ = x) solutions 
to Eq. (1.2). The periodic boundary conditions (2.14) 
translate into 



L 



26 (t) 



,v,f =/ 



28(t) 



v + a Le x ,t 



(2.15) 



In the case of the USF, these are the well-known Lees- 
Edwards's boundary conditions [71]. 

While the forms (1.2) and (2.10) of the Boltzmann 
equation arc fully equivalent, as are the respective bound- 
ary conditions (2.15) and (2.14), it is obvious that Eqs. 
(2.10) and (2.14) are much simpler to implement in com- 
puter simulations than Eqs. (1.2) and (2.15). This is 
especially important if one restricts oneself to uniform 
solutions of the form (2.13). In that case, Eq. (2.10) be- 
comes 



a/ 

dt 



The corresponding energy balance equation is 
2a 



Tit) 



(ITi 



P x i(t)-Cit)T(t), 



(2.17) 



where 



Note that 



dnT 



dvv 2 J[fJ}. 



(2.18) 



C(t) = 6(t)C(t) (2.19) 

and, thus, Eqs. (1.14) and (2.17) are equivalent. Al- 
though T = f, in Eq. (2.17) we keep the notation f to 
emphasize that here the temperature is seen as a function 
of the scaled time t. 

In cooling situations, i.e., in the USF if £ > 2\aP xy \/dp, 
in the ULF with a < if C > 2\a\P xx /dp, or in the ULF 
with ao > 0, the temperature can reach values much 
smaller than the initial one, which can cause technical 
difficulties (low signal-to-noise ratio) in the simulations. 
This is especially important in the ULF with ao > 
since no steady state exists and the temperature keeps 
decreasing without any lower bound. To manage this 
problem, it is convenient to introduce a velocity rescaling 
(or thermostat). From a mathematical point of view, let 
us perform the additional change of variables 



t = ti(t), 



7(v,?) = U(t)] fiv,t), 



(2.20) 



(2.21) 



where so far #(£) is an arbitrary (positive definite) proto- 
col function. The following identities are straightforward, 



-id/ 
dt 



2i 

dt 

d-l 



tf(t) d 



m j[fj} = j[fj]. 



(2.23) 



Therefore, Eq. (2.16) becomes 
df ^ d ~\ ^ d 

where 



a(t) ee 4^, 

m 



ft® 



(2.25) 



The Boltzmann equation [Eq. (2.24)] represents the 
action of a nonconservative external force F(v, t) = 
—ma(t)v(e x — m/x(t)v. The relationship between the 
granular temperatures defined from / and /, respectively, 

2 



is T(t) = T(t)/ 



1/2 



Thus, the thermostat choice 



keeps the rescaled temperature Tit) 

constant. While, at a theoretical level, Eq. (2.16) is sim- 
pler and more transparent than Eq. (2.24), the latter is 
more useful from a computational point of view in cooling 
situations. 
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B. Rheological functions 

In order to characterize the non-Newtonian properties 
of the unsteady USF and ULF, it is convenient to intro- 
duce generalized transport coefficients. 

As is well known, the NS shear viscosity 77NS is defined 
through the linear constitutive equation 



Pij = p5ij - ?/ns + Vj-Uj - -V ■ uSij j , (2.26) 

where we have taken into account that the NS bulk vis- 
cosity is zero in the low-density limit [14, 73]. This is 
especially relevant in compressible flows like the ULF. 
Making use of Eq. (1.1), we define (dimensionless) non- 
Newtonian viscosities rj*(a*) for the USF and the ULF 
by the relation 

p ij( a *) = % _ i]*(a*)a* y5 ix 5jg + 5j X Su - -5 x i5iA . 

(2.27) 

More specifically, setting ij — x£, Eq. (2.27) yields 

5 xt - P* e (a*) 



a* (1 + ^8 xt ) ' 



(2.28) 



The rheological function f]*(a*) differs in the USF from 
that in the ULF. In the latter flow it is related to the 
normal stress P xx . In that case, by symmetry, one has 

p* x + (d- i)p; y = d, so 

(2.29) 



xx yy 



-2r/*(a*)a* 



Since < P* x < d, the viscosity rj* (a*) in the ULF must 
be positive definite but upper bounded: 

7f(a * }< f^F, a*>0, (23[)) 

^2|a*|' a < U. 

In contrast, in the USF the viscosity function is related 
to the shear stress P* y , 



P 



-r]*(a*)a* 



(2.31) 



In this state the normal stress differences are character- 
ized by the viscomctric functions 



*5(o*) 



P* v {a*)-P* x {a*) 



7 *2 



PL^*)-P*{a*) 



(2.32) 
(2.33) 



III. KINETIC MODEL AND NONLINEAR 
HYDRODYNAMICS 

A. Kinetic model 



extension of the Bhatnagar-Gross-Krook (BGK) kinetic 
model [74]. in which the (inelastic) Boltzmann collision 
operator J[/, /] is replaced by a simpler form [75, 76]: 

J[f, f] -> -P{a)v{f - / hcs ) + £ A . [( V _ U )/]. (3.1) 

Here /hcs is the local version of the homogeneous cooling 
state distribution [12] and 



(d + 2)r(d/2)' 



d-X 



(3.2) 



is a convenient choice for the effective collision frequency. 
The factor /3(a) can be freely chosen to optimize agree- 
ment with the original Boltzmann equation. Although it 
is not necessary to fix it in the remainder of this section, 
we will take 



/3(a) = 1(1 + a) 



(3.3) 



at the end [76]. The cooling rate is defined by Eq. (1.11) 
but here we will take the expression obtained from the 
Maxwellian approximation, namely 



d + 2 
Ad 



(1 - a 2 )u. 



(3.4) 



This is sufficiently accurate from a practical point of view 
[55], especially at the level of the simple kinetic model 
(3.1). 

Using the replacement (3.1), Eq. (2.16) becomes 



Of d 
— - -a>o-=- 

dt dv x 



(vj) 



-M/-/hcs) + ^-(v/), (3.5) 



where v and C, are given by Eqs. (3.2) and (3.4), respec- 
tively, except for the change n — > n (recall that T = T). 
Taking second-order velocity moments on both sides of 
Eq. (3.5) one gets 



Pi 



-a \ PjiS ix + Ptf5 jx \ - C Pij ( 



p. 

1 1. 1 



nT8i 



(3.6) 

From the trace of both sides of Eq. (3.6) we recover 
the exact energy balance equation (2.17). The advan- 
tage of the BGK-like model kinetic equation (3.5) is that 
it allows one to complement Eq. (2.17) with a closed 
set of equations for the elements of the pressure tensor 
[17, 67]. It is interesting to note that Eq. (3.6), with 
/3 = (l + a)[d+ 1 + (d— l)a]/4d, can also be derived from 
the original Boltzmann equation in the Grad approxima- 
tion [19, 20, 53]. 

As discussed in Sec. I [cf. Eq. (1.19)], the relevant quan- 
tity is the reduced pressure tensor defined as 



In order to progress on the theoretical understanding 
of the USF and the ULF, it is convenient to adopt an 



Pj(t) p l3 {t) 



n{t)T(t) nT(t) 



(3.7) 
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Combining Eqs. (2.17) and (3.6) we obtain 



R 



■Pin 



where, according to Eq. (1.16), 



a *(J) = ^ = J?2_ 



(3.8) 



(3.9) 



is the reduced shear (£ = y) or longitudinal (£, = x) rate. 
Taking ij = x£ and ij = ££, Eq. (3.8) yields 



p* 

-^L = - a * (PI 
v 



2a* ? 

PM + — {P* t f 



(3.10) 



In contrast, the solution for the ULF (£ = x) is 

. = d(* fi + C 
Qs 2 f3 + d(*' 



p* 

xx.s 



f3 + dC 

P + C* 



(3.17) 



(3.18) 



Note that the steady-state values are independent of q. A 
linear stability analysis in the case of the USF [53] shows 
that the steady state, Eqs. (3.15) and (3.16), is indeed a 
stable solution of Eqs. (3.10)-(3.12). The proof can be 
easily extended to the ULF. 



C. Unsteady hydrodynamic solution 



?k = -2a*P* u 5 xl + 2 -?-P* u P* 



P(P* u -l). (3.11) 



Note that Eq. (3.11) is identical to Eq. (3.10) in the ULF 
(£ = x). Equations (3.10) and (3.11) must be comple- 
mented with the evolution equation for a*. We recall 
that v oc T q and a* oc T~ q with q = i. Thus, using Eq. 
(2.17) one simply obtains 



a 
v 



qa 



2a* 



where, according to Eq. (3.4), 
d + 2 



4d 



1 



a 



(3.12) 



(3.13) 



Here we will temporarily view q as a free parameter, so 
the solutions to Eqs. (3.10)-(3.12) depend parametrically 
on q. The exponent q is directly related to the wide class 
of dissipative gases introduced by Ernst el al. [77-79]. 

Equations (3.10)-(3.12) constitute a closed set of non- 
linear equations for {P* e (t), P^(€), a* (t)} that can be nu- 
merically solved subject to a given initial condition 



/o(v) => {PxiMi PlLQi a o}- 



Steady state 



(3.14) 



In the USF (£ = y), Eq. (3.8) implies that (P* z - 
P; y )/v = (2a*P* y /d - p)(P* zz - P; y ). Since, on physical 
grounds, a* P* y < 0, we conclude that P* z - P* y = in 
the hydrodynamic regime. Therefore, according to the 
kinetic model description, the second viscometric func- 
tion identically vanishes, i.e., 



'I' 



0. 



(3.19) 



Next, by symmetry, P* x + P* y + (d - 2)P* Z = d. This 
mathematical identity, combined with P zz = P yy , allows 
one to rewrite Eq. (2.32) as 



1 - P* 
** = -d M. 



(3.20) 



As sketched in Fig. 1 and described by Eqs. (1.17) 
and (1.19), the hydrodynamic solution requires the whole 
time dependence of P*j to be captured through a depen- 
dence on a* common to every initial state. As a first step 
to obtain such a hydrodynamic solution, let us eliminate 
time in favor of a* in Eqs. (3.10) and (3.11) with the help 
of Eq. (3.12), i.e., 



q[ i F * i + Q ) da* 



±00 



p:m + 2 (P* e y 



a 



8xi) ■ 



(3.21) 



Setting P* e = 0, P; e = 0, and a* = in Eqs. (3.10)- 
(3.12), they become a set of three (USF) or two (ULF) 
independent algebraic equations whose solution provides 
the steady-state values. In the case of the USF (£ = y) 
the solution is 



(3.15) 



2a* 

q[ —p* e + C 



OPI 

da* 



= —2Pl e 5 X £ 



-r nor „.o 



-£(P£-1). (3.22) 
a* 

This set of two nonlinear coupled differential equations 
must be solved in general with the initial conditions stem- 
ming from Eq. (3.14), namely 



„* ^ r p* p* p* p* i 

— "0 ^ \ r xt — r xl$i r U — r UfiS- 



(3.23) 



P* 

xy,s 



df3(* sgn(a:) 
2 P + C* ' 



P, 



ft 



P + C 



(3.16) 



Equations (3.21) and (3.22) must be solved in agreement 
with the physical direction of time. This means that 
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the solutions uncover the region |ajj| < |a*| < |a*| in 
conditions of cooling and the region |oq| > \a*\ > \a* s \ in 
conditions of heating (see Table I). In the case of the ULF 
with ciq > 0, one has Oq < a* < oo due to the absence of 
steady state. Equations (3.21) and (3.22) describe both 
the kinetic and hydrodynamic regimes. In order to isolate 
the hydrodynamic solution, one must apply appropriate 
boundary conditions [53, 67]. 

An alternative route to get the hydrodynamic solu- 
tion consists of expanding P* t and P£ e in powers of a* 
(Chapman— Enskog expansion), 

oc oo 

P* e (a*) =<y + £ c$a* j , P* u {a*) = 1 + £ cga**. 

3=1 i=i 

(3.24) 

Inserting the expansions in both sides of Eqs. (3.21) and 

(3.22) one can get the Chapman-Enskog coefficients c^J 

(i) 

and c ee in a recursive way. The first- and second-order 
coefficients are 



„(i) _ d+(d-2)5 xi _ (1) 

- d (/3 + qC) ' u - c *e° xi > 



(2) 

C xl 



c {2) 5 



(2 ) = 2(d - l)(d - 2 + g)6 xe + d(8 xl - 1) 
d 2 {P + qC)(/3 + 2qC) 



(3.25) 



(3.26) 



(3.27) 



Equation (3.25) gives NS coefficients, while Eqs. (3.26) 
and (3.27) correspond to Burnett coefficients. From Eqs. 
(3.24) and (3.25), Eq. (2.28) yields 



1 



lim ri* (a*) = 
a^O 1 v J /3 + q(* 



(3.28) 



Thus, the NS viscosity coincides in the USF and in the 
ULF, as expected. Regarding the USF first viscometric 
function, Eqs. (3.20), (3.24), and (3.27) gives the Burnett 
coefficient 



region. Second, even if \a*\ < |a*|, closed expressions for 
P* e {a*) and -P«(a*) are not possible. 

In order to get closed and explicit (albeit approximate) 
solutions, we formally take q as a small parameter and 
perturb around q = [62, 66], 



P* j (a*) = P*V(a*)+qP*f L) (a*) + ---- 
Setting q = in Eqs. (3.21) and (3.22) one gets 



p:i 0) = d 



"(0) 



1 



1 + 2j e (a*) ' 



(3.30) 



(3.31) 



(3.32) 



where 7^(a*) is the physical solution of a cubic (USF, 
I = y) or a quadratic (ULF, I — x) equation: 



The respective solutions are 
1 



ly{a*) = ^sinh 2 



7y (l + 2 7y ) 2 = — , 



d(l-%](l + 27x ) = l. 



*2 



cosh" 1 f 1 + 27° 



2 d 



(3.33) 



(3.34) 



, (3.35) 



, a* 1 1 fa* 1\ 2a* 

lx{a)= 2p~A + 2\\J + 2) ~W 



For small \a*\ one has 



,*2 



7«(o*) = ^ 5 + 0(a* a ), 



(3.36) 



(3.37) 



lim \l>T(a*) 

a*— >0 



(P + q(*)(l3 + 2q(*)' 



(3.29) 



In general, all the even (odd) coefficients of P* y (Py V ) 
vanish in the USF. In the ULF, however, all the coef- 
ficients of P* x are non-zero. It is interesting to remark 
that, in contrast to the elastic case (£* = 0) [66, 80], 
the Chapman-Enskog expansions (3.24) are convergent 
[59, 60, 67] if > 0. On the other hand, the radius 
of convergence is finite and coincides with the stationary 
value |a*|. 

The series (3.24) clearly correspond to the hydrody- 
namic solution since they give P* e and P£ e as unambigu- 
ous functions of a* , regardless of the details of the initial 
conditions (3.23). However, the series have two short- 
comings. First, since they diverge for |o*| > |a*|, they do 
not provide P* t (a*) and P^ t {a*) in a direct way for that 



so 

D *(o)/ *\ c d + (d — 2)8 x i * ~i ( *2n Qft i 
P xt \ a ) = °ni a )i (3.39) 



in agreement with Eq. (3.25). Furthermore, it is easy to 
check that in the steady state [cf. Eqs. (3.15) and (3.17)] 
one gets 



l c* 

7v(*.) = 2j. 



lx{a* s ) 



d-1 C 
2~ [3 + dC ' 



(3.40) 



(3.41) 
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so Eqs. (3.31) and (3.32) yield Eqs. (3.16) and (3.18), as 
expected. 

Once P*^ are known, Eqs. (3.21) and (3.22) provide 
P*j ■ In the USF case (£ = y) the results arc 

= ~P^ ] h y {a*\ (3.42) 



p;« = 6P<°) ly (a*)H y (a*), (3.43) 



where 



{al s fCV/?-^(g*)][i- 6 7 >*)] 

[l + 6 7 ,(a*)] 2 



^)^ W - 27> ; ) , (3-45) 
[l + 6 7 ,(a*)] 2 

and use has been made of Eq. (3.33) and the relation 
dl v {a*) 2a* 1 



da* p 2 d[l + 2 ly {a*)][l+Q ly {a*)\ 



(3.46) 



Analogously, taking into account Eq. (3.34), the final ex- 
pression in the ULF case {I = x) is 



P:W=2P<°h x (a*)h x (a*), 



(3.47) 



where 



M**)= 2aV ^ 27 * (a * ) + C T (3-48) 
[l-2a*/P + 4 lx (a*)} 2 

Note that in the steady state P^ = both in the 
USF and the ULF. This is consistent with the fact that 
the steady-state values are independent of q. For small 
lo*|. 



* N _ d+(d-2)S x£ / 
p 2 d 



^TV) = - ' ^f^ Ca* + 0(an, (3-49) 



which again agrees with Eq. (3.25). 

In principle, it is possible to proceed further and get 
the terms of orders q 2 , q 3 , . . . , in Eq. (3.30) [66]. How- 
ever, for our purposes it is sufficient to retain the linear 
terms only. 



D. Rheological model 

The definitions of the rheological functions (2.28), 
(2.32), and (2.33) are independent of any specific model 
employed to obtain P*j(a*). Here we make use of the 
kinetic model (3.1) and the expansion (3.30) truncated 
to first order in q. Next, by means of a Pade approxi- 
mant we construct (approximate) explicit expressions for 



rj*(a*). Let us start with the ULF (i = x), in which case 
Eqs. (3.31) and (3.47) yield 



rj*(a*) 



lx(a*) 



1 



d- 1 a*[l + 2 lx {a*)] 1 + qh x {a*) ' 



(3.50) 



Analogously, in the USF case 
and (3.42) give 

V *(a*) 



p[l + 2 ly (a*)f l + qh v {a*Y 



y), Eqs. (3.31), (3.33), 



(3.51) 



Let us analyze now the USF first viscometric function. 
Using Eqs. (3.32), (3.33), and (3.43), Eq. (3.20) gives 



/3 2 [1 + 2 7y (a*)r 



[l-3qH y (a*)]+0(q 2 ). 

(3.52) 

Again, it is convenient to construct a Pade approximant 
of **(a*). Here we take 



^{(a*) ps -- 



1 



P 2 [1 + 2 7y (a*)] 3 [1 + qH y {a*)][l + 2qH y (a*)} 

(3.53) 

In principle, one should have written 1 + 3qH y instead of 
(1 + qH y )(l + 2qH y ) in Eq. (3.53), but the form chosen 
has the advantage of being consistent with the Burnett 
coefficient (3.29) for any q. 

In summary, our simplified rheological model consists 
of Eq. (3.50), complemented with Eqs. (3.36) and (3.48), 
for the ULF and Eqs. (3.51) and (3.53), complemented 
with Eqs. (3.35), (3.44), and (3.45), for the USF. Since 
we are interested in hard spheres, we must take q = i in 
those equations. 

This approximation has a number of important prop- 
erties. First, as said in connection with the Chapman- 
Enskog expansion (3.24), Eqs. (3.50), (3.51), and (3.53) 
qualify as a (non-Newtonian) hydrodynamic description. 
Second, in contrast to the full expansions (3.24) and 
(3.30), they provide the relevant elements P* g of the pres- 
sure tensor as explicit functions of both the shear or lon- 
gitudinal rate a* and the coefficient of normal restitution 
a (through j3 and (*). Third, as seen from Eq. (3.28), 
Eqs. (3.50) and (3.51) agree with the exact NS coeffi- 
cients predicted by the kinetic model (3.1) for arbitrary 
values of the parameter q; this agreement extends to the 
Burnett-order coefficient (3.29). Next, the correct steady- 
state values [Eqs. (3.15)-(3.18)] are included in the Pade 
approximants (3.50), (3.51), and (3.53). Finally, it can 
be checked that the correct asymptotic forms in the limit 
|a*| — > oo [53, 67] are preserved. 

Moreover, even in the physical case of hard spheres 
(q = ^), Eqs. (3.50), (3.51), and (3.53) represent an ex- 
cellent analytical approximation to the numerical solu- 
tion of the set of Eqs. (3.21) and (3.22) with appropriate 
initial conditions [53, 67]. The results for a = 0.5 and 
a = 0.9 are presented in Figs. 4 and 5 in the cases of the 
USF and the ULF, respectively. We observe a generally 
good agreement between the numerical and the simpli- 
fied results. This is especially true in the USF, where the 
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FIG. 4. (Color online) Shear-rate dependence of (a) the 
viscosity function and (b) the first viscometric function for 
a = 0.5 and a = 0.9 in the USF with d = 3. The solid lines 
have been obtained from numerical solutions of Eqs. (3.21) 
and (3.22), while the dashed lines correspond to the simpli- 
fied rheological model (3.51) and (3.53). The circles represent 
the steady-state values [cf. Eqs. (3.15) and (3.16)]. Note that 
the simplified solution deviates practically from the numerical 
solution only for the most inelastic system (a = 0.5) and in 
the region a* < a*. 



limitations of the rheological model are only apparent for 
the most inelastic system (a = 0.5) and for shear rates 
smaller than the steady-state value, the agreement be- 
ing better for the viscosity than for the first viscometric 
function. In the ULF case the differences are more im- 
portant, both for a = 0.5 and a = 0.9, although they are 
restricted to longitudinal rates near the maximum and 
also to values more negative than that of the maximum. 



IV. SIMULATION DETAILS 

We have performed DSMC simulations of the three- 
dimensional (d = 3) USF and ULF. The details are sim- 
ilar to those described elsewhere [55] , so here we provide 
only some distinctive features. 

Since the original Boltzmann equation (1.2) is fully 
equivalent to the scaled form (2.10) [with the change of 
variables (2.1)-(2.5)], we have solved the latter and have 
applied the periodic boundary conditions (2.14). We call 
this the "inhomogeneous" problem since the scaled dis- 
tribution function f{£, v, t) is, in principle, allowed to de- 
pend on the scaled spatial variable £. On the other hand, 
if one restricts oneself to uniform solutions (2.13), then 
Eq. (2.10) reduces to Eq. (2.16). The solution to Eq. 
(2.16) by the DSMC method will be referred to as the 
"homogeneous" problem. Most of the results that will 



1.5 - 
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0.5 - 
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FIG. 5. (Color online) Longitudinal-rate dependence of the 
viscosity for a = 0.5 and a = 0.9 in the ULF with d — 3. The 
solid lines have been obtained from numerical solutions of Eq. 
(3.21), while the dashed lines correspond to the simplified 
rheological model (3.50). The circles represent the steady- 
state values [cf. Eqs. (3.17) and (3.18)]. 



be presented in Sec. V correspond to the homogeneous 
problem. 

A wide sample of initial conditions fo{£,v) has been 
considered, as described below. Note that, since 9(0) = 1 
[cf. Eq. (2.3)], 1 = £ and v = v - a £e x at t = 0. 
Consequently, n (£) = n (£) and u (£) = u (£) — ao£e x . 
However, for consistency, we will keep the tildes in the 
expressions of the initial state. An exception will be the 
initial temperature because T = T for all times. 

The inhomogeneous problem for the USF was already 
analyzed in Ref. [55] and, thus, only the inhomogeneous 
problem for the ULF is considered in this paper. The 
chosen initial condition is 

/ (x,v) = ^^(|v-Uo(^|-v^oM)> (4-1) 

where 5(x) is Dirac's distribution and the initial density 
and velocity fields are 



n (x) 



Uq(x) = CLqL cos 



1 2ttx\ 
l + -sm— j, 



L 



(4.2) 



(4-3) 



respectively, while the initial temperature To is uniform. 
In Eq. (4.2) (n) is the density spatially averaged between 
x = —L/2 and x = L/2. This quantity is independent of 
time. In our simulations of the inhomogeneous problem 
we have taken a = 0.5 for the coefficient of restitution, 



11 



TABLE II. Values of the initial pressure tensor for the four 
initial conditions of the class (4.6). 



Label 





Pxx,0 


Pyy>o 


Pzz,0 




BO 





2nT 





nT 





Bl 


tt/4 


nT 


nT 


nT 


— nT 


B2 


7r/2 





2nT 


nT 





B3 


3tt/4 


nT 


nT 


nT 


nT 



ao = — 4/to for the initial longitudinal rate, and L 
for the (scaled) box length. Here, 



A = 



1 



A 



V2ir(n)a 2 



to = 



2.5A 



(4.4) 



are a characteristic mean free path and an initial charac- 
teristic collision time, respectively. 

As for the homogeneous problem (both for USF and 
UFL) two classes of initial conditions have been chosen. 
First, we have taken the local equilibrium state (initial 
condition A), namely 



/o(v) = n 



( m 



3/2 



2 /2T 



(4.5) 



The other class of initial conditions is of the anisotropic 
form 

1/2 



A(v) 



(-. 



Ill 



n 

2 V27rTo 

x [S (v x - 
+5 (v x 



K/2T 

- Vq cos 4») 5 (vy - 
Vq cos </>) 5 (v y — 



- Vq sin 4>) 
V sin 0)], (4.6) 



where Vq = 2Tq /m is the initial thermal speed and <fi € 
[0,7r] is an angle characterizing each specific condition. 
The pressure tensor corresponding to Eq. (4.6) is given 
by Pxxfl = 2nT cos 2 (f>, P yyfl = 2nT sin 2 (f>, P zzfi = nT , 
and P X yfl = —nTo sin 2(f). The four values of </> considered 
are <f> = kir/4 with k = 0, 1, 2, and 3; we will denote the 
respective initial conditions of type (4.6) as BO, Bl, B2, 
and B3. The values of the elements of the pressure tensor 
for these four initial conditions are displayed in Table II. 
In the fully developed USF it is expected that P xy < 
(if a > 0) and P xx > P yy . As we see from Table II, 
the four initial conditions are against those inequalities, 
especially in the case of condition B3. As for the fully 
developed ULF, the physical expectations are P xy = 0, 
Pyy = Pzz, and P xx < P yy if a > and P xx > P yy if 
ao < [cf. Eq. (2.29)]. Again, none of the four initial 
conditions is consistent with those physical expectations, 
especially in the case of condition BO if ao > and B2 if 
ao < 0. The "artificial" character of the initial conditions 
(4.6) represents a stringent test of the scenario depicted 
in Fig. 1. 

As summarized in Table I, when £ > 2\aoP xy \/3nT in 
the USF or when C > -2a P xx /3nf in the ULF, the tem- 
perature decreases with time (cooling states) either with- 
out lower bound (ULF with ao > 0) or until reaching the 



steady state (ULF with ao < and USF). In those cool- 
ing states the temperature can decrease so much (relative 
to the initial value) that this might create technical prob- 
lems (low signal-to-noise ratio), as mentioned at the end 
of Sec. II. This can be corrected by the application of a 
thermostatting mechanism, as represented by the change 

of variables (2.20) and (2.21) with tf(t) cx T(t) 
The DSMC implementation of Eqs. (2.24) and (2.25) is 
quite simple. Let us denote by {vj(t); i = 1, . . . , N} the 
(rescaled) velocities of the N simulated particles at time 
t. The corresponding rescaled temperature and shear (or 
longitudinal) rate are T(t) and a(t), respectively. Dur- 
ing the time step St the velocities change due to the ac- 
tion of the (deterministic) nonconservative external force 
—ma(t)v£e x and also due to the (stochastic) binary col- 
lisions. Let us denote by {v£(f + St); i = 1, . . . , N} and 
by T'(t + St) the velocities and temperature after this 
stage. Thus, the action of the thermostat force —mfx(t)v 
is equivalent to the velocity rescaling v£(i + St) — > 

9,(?+ St) = v[(t + St)^f(t)/f'(t + St), so f'{t + St) -> 
T(t + St) = T(t). Similarly, the rescaled shear or longitu- 



dinal rate is updated as a(t+St) = a(t)yT'(t + St)/T(t), 

so a (t + St) I ^Jf'it + St) = a (t) / \Jf^). 

For each one of the five initial conditions for the homo- 
geneous problem we have considered three coefficients of 
restitution: a = 0.5, 0.7, and 0.9. In the case of the USF, 
the values taken for the shear rate have been a = 0.01/t , 
a = 0.1/to, a = 4/to, and a = 10/ro. The two first val- 
ues (a = 0.01/to and a = 0.1/ro) are small enough to 
correspond to cooling cases, even for the least inelastic 
system (a = 0.9), while the other two values (a = 4/to 
and a = 10/to) are large enough to correspond to heat- 
ing cases, even for the most inelastic system (a = 0.5). 
In the case of the ULF we have chosen ao = 0.01/ro, 
ao = —0.01/ro, and ao = —10/ro. The first and second 
values correspond to cooling states (without and with a 
steady state, respectively), while the third value corre- 
sponds to heating states. Therefore, the total number 
of independent systems simulated in the homogeneous 
problem is 60 for the USF and 45 for the ULF. 

The technical parameters of the simulations have been 
the following ones: TV = 10 6 simulated particles, an adap- 
tive time step St = 10~ 3 to\/Tq/ (T) , and a layer thick- 
ness (inhomogeneous problem) Sx = 0.05A. Moreover, 
in order to improve the statistics, the results have been 
averaged over 100 independent realizations. 



V. RESULTS 

A. USF. Homogeneous problem 

We have simulated the Boltzmann equation describing 
the homogeneous problem of the USF, i.e., Eq. (2.16) 
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with £ = y, by means of the DSMC method. As said in 
Sec. IV, three coefficients of restitution (a = 0.5, 0.7, and 
0.9) and four shear rates (a = 0.01/ro, 0.1/to, 4/to, and 
10/to) have been considered. For each combination of a 
and a, five different initial conditions (A and B0-B3) have 
been chosen. In the course of the simulations we focus 
on the temporal evolution of the elements of the reduced 
pressure tensor P*j, Eq. (3.7), and of the reduced shear 
rate a* , Eq. (3.9). From these quantities one can evaluate 
the viscosity 77*, Eq. (2.31), the first viscometric function 
"J*, Eq. (2.32), and the second viscometric function 
Eq. (2.33). The effective collision frequency v in Eq. (3.9) 
is defined by 
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The factor 1.016 comes from an elaborate Sonine approx- 
imation employed to determine the NS shear viscosity 
t/ns °f a S as °f elastic hard spheres [14]. For simplicity, 
and to be consistent with the approximate character of 
the kinetic model (3.1), this factor is not included in Eq. 
(3.2). 

Note that the initial reduced shear rate is Oq = 
aro/0.8885. Time is monitored through the accumulated 
number of collisions per particle, i.e., the total number 
of collisions in the system since the initial state, divided 
by the total number of particles. The reduced quantities 
and the number of collisions per particle are not affected 
by the changes of variables discussed in Sec. II. 

As a representative case, we first present results for the 
most inelastic system (a = 0.5). Figures 6 and 7 show 
the evolution of a* and 77* for the cooling states (a = 
0.01/ro and 0.1/to) and the heating states (a = 4/to and 
10/to), respectively. We clearly observe that after about 
30 collisions per particle (cooling states) or 20 collisions 
per particle (heating states) both a* and 77* have reached 
their stationary values. 

Figure 6 shows that, for each value of a, the full tem- 
poral evolution of a* cx T^ 1 / 2 is practically independent 
of the initial condition, especially in the case a = 0.01/ro. 
This is due to the fact that for these low values of ar 
the viscous heating term —2aP xy /dn in Eq. (1.14) can be 
neglected versus the inelastic cooling term (T for short 
times, so the temperature initially evolves as in the ho- 
mogeneous cooling state (decaying practically exponen- 
tially with the number of collisions), hardly affected by 
the details of the initial state. On the other hand, the 
first stage in the evolution of the reduced viscosity 77* is 
widely dependent on the type of initial condition, as ex- 
pected from the values of P xy ,o shown in Table II. In the 
heating cases, Fig. 7 shows that the evolution followed 
by both a* and 77* is distinct for each initial condition, 
except when the steady state is practically reached. 

In any case, the interesting point is whether an un- 
steady hydrodynamic regime is established prior to the 
steady state. If so, a parametric plot of 77* versus a* 
must approach a well-defined function 77* (a*), regardless 
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FIG. 6. (Color online) (a) Reduced shear rate a* and (b) re- 
duced viscosity 77* versus the number of collisions per particle 
for the USF with a = 0.5 in the cooling states a = 0.01/ro 
[blue (dark gray) lines] and a = 0.1/ro [orange (light gray) 
lines]. The legend refers to the five initial conditions consid- 
ered. The dotted horizontal line in panel (a) denotes the value 
a£ = 0.4 above which the hydrodynamic regime is clearly es- 
tablished (see text). 
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FIG. 7. (Color online) (a) Reduced shear rate a* and (b) re- 
duced viscosity 77* versus the number of collisions per particle 
for the USF with a — 0.5 in the heating states a — 4/ to [or- 
ange (light gray) lines] and a — 10/rn [blue (dark gray) lines]. 
The legend refers to the five initial conditions considered. The 
dotted horizontal line in panel (a) denotes the value = 1.25 
below which the hydrodynamic regime is clearly established 
(see text). 
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FIG. 8. (Color online) (a) Reduced viscosity rf , (b) first 
viscometric function — ty*, and (c) second viscometric func- 
tion ^2 versus the reduced shear rate a* for the USF with 
a — 0.5 in the cooling states a = 0.01/ro [blue (dark gray) 
lines] and a = 0.1/to [orange (light gray) lines]. The legend 
refers to the five initial conditions considered. The circles 
represent the steady-state points (a*, 77*), (aj,— ^i^), and 
(ets,^^), respectively. The dotted vertical line denotes the 
value a\ = 0.4 above which the curves collapse to a common 
one. 



FIG. 9. (Color online) (a) Reduced viscosity rf , (b) first 
viscometric function — and (c) second viscometric func- 
tion ^2 versus the reduced shear rate a* for the USF with 
a — 0.5 in the heating states a — 10/ro [blue (dark gray) 
lines] and a — 4/r [orange (light gray) lines]. The legend 
refers to the five initial conditions considered. The circles 
represent the steady-state points (a*,rj*), (o*,-$J ,), and 
(aj,$2,s); respectively. The dotted vertical line denotes the 
value a£ = 1.25 below which the curves collapse to a common 
one. 



of the initial condition. Figures 8 and 9 present such a 
parametric plot, also for the viscometric functions, for 
the cooling and heating states, respectively. We observe 
that, for each class of states (either cooling or heating) 
, the 10 curves arc attracted to a common smooth "uni- 
versal" curve, once the kinetic stage (characterized by 
strong variations, especially in the case of the second vis- 
cometric function for the cooling states) is over. One can 
safely say that the hydrodynamic regime extends to the 
range 0.4 < a* < a* for the considered cooling states and 
to the range 1.25 > a* > a* for the considered heating 
states. We will denote the above threshold values of a* 
by a£. It is expected that a£ depends on the initial value 
dp (apart from a weaker dependence on the details of the 
initial distribution); in fact, Figs. 8 and 9 show that the 
hydrodynamic regime is reached at a value a* h < 0.4 by 
the states with a = 0.01/ro and at a value a* h > 1.25 
by the states with a = 10/ro- Here, however, we adopt 
a rather conservative criterion and take a common value 
a* h = 0.4 for a — 0.01/ro and O.lro and a common value 



a* h = 1.25 for a = 4/to and 10to. It is also quite appar- 
ent from Figs. 8 and 9 that the collapse to a common 
curve takes place earlier for rj* than for ^l, ^2 being the 
quantity with the largest "aging" period. 

From Fig. 6 it can be seen that the value = 0.4 is 
reached after about 5 collisions per particle in the states 
with a = 0.1/to and after about 15 collisions per par- 
ticle in the states with a = 0.01/to. Similarly, Fig. 7 
shows that the value a£ = 1.25 is reached after about 
5 collisions per particle in the states with a = 4/ro and 
a = 10/ro. Given that, as said before, the values a ; * = 0.4 
and a* h = 1.25 are conservative estimates, we find that, 
as expected, the duration of the kinetic stage is shorter 
than the duration of the subsequent hydrodynamic stage, 
before the steady state is eventually attained. 

We have observed behaviors similar to those of Figs. 
6-9 for the other two coefficients of restitution (a = 0.7 
and a = 0.9, not shown). Table III gives the values of 
and the duration of the aging and transient periods 
for the 12 classes of states analyzed. It turns out that 
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TABLE III. This table shows, for each value of the coefficient 
of restitution a and each value of are, the duration of the ag- 
ing period toward the unsteady hydrodynamic regime and the 
total duration of the transient period toward the steady state, 
both measured by the number of collisions per particle. The 
aging period is defined as the number of collisions per particle 
needed to reach a common threshold value a* h above (below) 
which the hydrodynamic regime is established for the cooling 
(heating) states. The table also includes the stationary value 
a* of the reduced shear rate. 
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the number of collisions per particle the system needs to 
lose memory of its initial state is practically independent 
of a for the heating states. However, the total duration 
of the transient period increases with a [58]. In fact, 
there is no true steady state in the elastic limit a — > 
1. Therefore, the less inelastic the system, the smaller 
the fraction of the transient period (as measured by the 
number of collisions per particle) spent by the heating 
states in the kinetic regime. In the cooling cases, both 
the aging and the transient periods increase with a. 

Figure 10 displays the viscosity 77* (a*) and the visco- 
metric functions -**(a*) and ^(a*) for a = 0.5, 0.7, 
and 0.9, both for the cooling and the heating states. Here 
we have focused on the ranges of a* where the hydro- 
dynamic regime can safely be assumed to hold, namely 
0.4 < a* < 1.25 for a = 0.5, 0.35 < a* < 1 for a = 0.7, 
and 0.2 < a* < 0.8 for a = 0.9. The curves describ- 
ing the predictions of the simplified rheological model for 
77*, Eq. (3.51), and for Eq. (3.53), are also included. 
We can see that the curves corresponding to the cooling 
states (a* < a*) and those corresponding to the heat- 
ing states (a* > a*) smoothly match at the steady-state 
point. In the case of the nonlinear viscosity function 77*, 
the 10 curves building each branch (cooling or heating) 
for each value of a exhibit a very high degree of over- 
lapping. Due to fluctuations associated with the normal 
stress differences, the common hydrodynamic curves for 
the viscometric functions are much more coarse grained, 
especially in the case of VP^i whose magnitude is at least 
10 times smaller than that of ^l- The impact of fluctu- 
ations is higher in the cooling branches (a* < a* s ) than 
in the heating branches (a* > a*). In fact, the definition 
of the viscometric functions [sec Eqs. (2.32) and (2.33)] 
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FIG. 10. (Color online) (a) Reduced viscosity 77*, (b) first 
viscometric function — and (c) second viscometric func- 
tion ^2 versus the reduced shear rate a* for the USF with, 
from top to bottom, a = 0.5, a = 0.7, and a = 0.9. In 
order to focus on the hydrodynamic regime, the curves have 
been truncated at the values of a* h given by Table III. The 
circles represent the steady-state points. The thin solid lines 
in panels (a) and (b) represent the predictions of our sim- 
plified rheological model, Eqs. (3.51) and (3.53). Note that 
the model is unable to predict a non-zero second viscometric 
function. 



shows that the signal-to-noise ratio is expected to dete- 
riorate as the reduced shear rate a* decreases. Since 77*, 
— and ^2 decrease with decreasing inelasticity, the 
role played by fluctuations increase as a increases. It is 
also interesting to remark that, despite its simplicity and 
analytical character, the rheological model described by 
Eqs. (3.51) and (3.53) describes very well the nonlinear 
dependence of 77* (a*) and ^l(a)*. On the other hand, 
the simple kinetic model (3.1) does not capture any dif- 
ference between the normal stresses P yy and P zz in the 
hydrodynamic regime and, thus, it predicts a vanishing 
second viscometric function [cf. Eq. (3.19)]. 

Figure 10 strongly supports Eq. (1.19) in the USF, i.e., 
the existence of well-defined hydrodynamic rheological 
functions P*j(a*) [or, equivalently, 77* (a*) and ^t t2 ( a *)] 
acting as "attractors" in the evolution of the pressure 
tensor Pjj(t|/o), regardless of the initial preparation /o. 
The stronger statement (1.17) (see Fig. 1) is also sup- 
ported by the simulation results [58]. 
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FIG. 11. (Color online) Profiles of (a) density, (b) flow ve- 
locity, and (c) temperature in the ULF, starting from the 
initial condition of Eqs. (4.1)-(4.3). The curves correspond 
to (1) t = 0.046-ro 0.03 coll/part), (2) t = 0.158r 
(» 0.29 coll/part), (3) t = 0.306r (w 0.72 coll/part), and 
(4) t = 0.400-to (» 1.07 coll/part). The dashed lines repre- 
sent the initial profiles. 



B. ULF. Inhomogeneous problem 

Now we turn to the ULF sketched by Fig. 3. By per- 
forming the changes of variables (2.1), (2.4), and (2.5) 
(with I = x), we have numerically solved the Boltzmann 
equation (2.10) by means of the DSMC method. Al- 
though Eq. (2.10) admits homogeneous solutions in the 
fully developed ULF [see Eq. (2.13)], it is worth checking 
that Eq. (2.10), complemented by the periodic boundary 
conditions (2.14), indeed leads an inhomogeneous initial 
state toward a (time-dependent) homogeneous state. A 
similar test was carried out in the case of the USF in Ref. 
[55]. 

As described in Sec. IV, we have considered the highly 
inhomogeneous initial state given by Eqs. (4.1)-(4.3) 
with a,Q = — 4/to and L = 2.5A, and solved Eq. (2.10) for 
a coefficient of restitution a = 0.5. The instantaneous 
density, flow velocity, and temperature profiles are plot- 
ted in Fig. 11 at four representative times. In order to 
decouple the relaxation to a homogeneous state from the 
increase of the global temperature (here viscous heating 
prevails over inelastic cooling), panel (c) of Fig. 11 dis- 



plays the ratio T/(T), with (T) = (p)/(n), where (n) and 
(p) are the density and hydrostatic pressure, respectively, 
averaged across the system. 

We observe that at t = 0.046ro the density, velocity, 
and temperature profiles arc still reminiscent of the initial 
fields, except in the region —0.4 <x/L< —0.3, where the 
density and the temperature present a maximum and the 
flow velocity exhibits a local minimum. At t = 0.158ro 
the inhomogencitics arc still quite strong, with a widely 
depopulated region 0.2 < x/L < 0.5 of particles practi- 
cally moving with the local mean velocity (which implies 
an almost zero temperature). By t = 0.306ro the pro- 
files have smoothed out significantly. Finally, the system 
becomes practically homogeneous at t = 0.400to. Thus, 
the relaxation to the homogeneous state lasts about 1 
collision per particle only. The system keeps evolving 
to the steady state, which requires about 20 collisions 
per particle. In fact, (T) ~ 58To after 1 collision per 
particle, while (T) ~ 211To in the steady state. Recall 
that 5(2;, t) = translates into u(ar, t) = a(t)xe x [see Eq. 
(2.12)]. _ 

It is important to remark that the relaxation to the 
ULF geometry observed in Fig. 11 does not discard the 
possibility of spatial instabilities for sufficiently large val- 
ues of the system size L, analogously to what happens in 
the USF case [27, 29, 31, 37-39, 81]. 



C. ULF. Homogeneous problem 

Now we restrict to the homogeneous ULF problem. 
The homogeneous Boltzmann equation (2.16) wit i = x 
has been solved via the DSMC method outlined in Sec. 
IV for a = 0.5, 0.7, and 0.9, starting from the initial con- 
ditions (4.5) and (4.6) with ao = — 10/ro, — 0.01/ro, and 
0.01/ro- Note that the Bl initial state becomes B3, and 
vice versa, under the change v y — ¥ ~v y , and so both are 
formally equivalent in the ULF geometry. As said before, 
a steady state is only possible if ao < 0. Moreover, the 
choices ao = ±0.01/to correspond to cooling states, while 
the choice ao = — 10/to corresponds to heating states. In 
the course of the simulations the reduced longitudinal 
rate a*, Eq. (3.9), and the reduced nonlinear viscosity 
rf, Eq. (2.29), are evaluated. 

As in the USF case, let us adopt a = 0.5 to illustrate 
the behaviors observed. Figure 12 displays the time evo- 
lution of | a* | and rj* for the 15 simulated states (5 for each 
value of ao). The states with negative ao become station- 
ary after about 20 collisions per particle, a value compa- 
rable to what we observed in the USF case for a = 0.5 
(see Table III). As for the states with ao = 0.01/ro, a* 
monotonically increases and rj* monotonically decreases 
(except for a possible transient stage) without bound. It 
is worth noticing that the first stage of evolution (up to 
about 7 collisions per particle) of |o* | and rj* for the initial 
condition B0 (B2) with ao = — 0.01/to is very similar to 
those for the initial condition B2 (B0) with ao = 0.01/tq. 
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FIG. 12. (Color online) (a) Absolute value of the reduced 
longitudinal rate \a*\ and (b) reduced viscosity r/* versus the 
number of collisions per particle for the ULF with a — 0.5 
in the cooling states eto = — 0.01/to [blue (dark gray) lines] 
and ao = 0.01/ro (black lines) and in the heating states ao = 
— 10/ To [orange (light gray) lines]. The legend refers to the 
five initial conditions considered. The dotted horizontal lines 
in panel (a) denote the values — ±0.08 and a£ = —0.4 (see 
text). 
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FIG. 13. (Color online) Reduced viscosity 77* versus the ab- 
solute value of the reduced longitudinal rate \a*\ for the ULF 
with a — 0.5 in the cooling states ao = —0.01/ro [blue (dark 
gray) lines] and ao = 0.01/ro (black lines) and in the heating 
states ao = — 10/ro [orange (dark gray) lines]. The legend 
refers to the five initial conditions considered. The circle rep- 
resents the steady-state point (|a*|,??*). The dotted vertical 
lines denote the values a^ = ±0.08 and a£ = —0.4 (see text). 



Eliminating time between a* and rj* one obtains the 
parametric plot shown in Fig. 13. In the cooling states 
ao = ±0.01/to we observe that the hydrodynamic regime 
is reached at approximately |a£| = 0.08. From Fig. 12 we 
see that this corresponds to about seven to eight collisions 
per particle. The heating states ao = — 10/tq deserve 
some extra comments. In those cases the time evolution 
is so rapid that, strictly speaking, the collapse of the five 
curves takes place only for a* > a* > a£ = —0.4, which 
corresponds to an aging period of about three to four 
collisions per particle (see Fig. 12). On the other hand, 
it can be clearly seen from Fig. 13 that the three curves 
corresponding to the initial states B0, Bl, and B3 have 
collapsed much earlier and are not distinguishable on the 
scale of the figure. It seems that the isotropic local equi- 
librium initial state A and the highly artificial anisotropic 
initial state B2 require a longer kinetic stage than in the 
cases of the initial states B0, Bl, and B3. While the rel- 
atively slower convergence of the initial condition B2 can 
be expected because of its associated "incorrect" negative 
viscosity, it seems paradoxical that the local equilibrium 
initial condition A also relaxes more slowly than the ini- 
tial conditions Bl and B3, the three of them having a 
zero initial viscosity. This might be due to the isotropic 
character of the local equilibrium distribution, in contrast 
to the high anisotropy of conditions Bl and B3. In what 
follows we will discard the initial conditions A and B2 for 
ao = — 10/to and assume that the states starting from the 
initial conditions BO, Bl, and B3 have already reached 
the hydrodynamic stage for say a* > —2. A stricter limi- 
tation to a* > —0.4 would miss the interesting maximum 
of rf at a* < a* predicted by the BGK-likc kinetic model 
(see Fig. 5). In any case, as observed in Fig. 13, the be- 
havior of the curves with ao = — 10/to corresponding to 
the initial conditions B2 and, especially, A is very close 
to the one corresponding to the initial conditions B0, Bl, 
and B3. 

The analysis for the cases a = 0.7 and 0.9 is similar 
to the one for a = 0.5 and, thus, it is omitted here. As 
in the USF (see Table III), the main effect of increasing 
a is to slow down the dynamics: the steady state (if 
a < 0) is reached after a larger number of collisions and 
the hydrodynamic stage requires a longer period. 

The a*-dependence of 77* for the three values of a is 
shown in Fig. 14, where we have focused on the intervals 
— 2 < a* < a* for a = — 10/r (initial conditions B0, Bl, 
and B3), -0.08 > a* > a* for a = -0.01/t (initial con- 
ditions A and B0-B3), and a* > 0.08 for a = 0.01/r 
(initial conditions A and B0-B3). Analogously to the 
case of Fig. 10, it can be seen that the heating and cool- 
ing branches with negative ao smoothly match at the 
steady state. Moreover, the cooling branch with posi- 
tive ao is a natural continuation of the cooling branch 
with negative ao, even though the zero longitudinal rate 
a* = represents a repeller in the time evolution of both 
branches. Figure 14 also includes the predictions of the 
rheological model (3.50). The agreement with the sim- 
ulation results is quite satisfactory, although the model 
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FIG. 14. (Color online) Reduced viscosity n* versus the re- 
duced longitudinal rate a* for the ULF with, from top to 
bottom, q = 0.5, a = 0.7, and a = 0.9. The circles represent 
the steady-state points. The thin dashed lines represent the 
predictions of our simplified rheological model, Eq. (3.50). 



tends to underestimate the maxima. This discrepancy is 
largely corrected by the true numerical solution of Eq. 
(3.21) (see Fig. 5). However, as done in Fig. 10, we pre- 
fer to keep the rheological model due to its explicit and 
analytical character. 



VI. CONCLUSIONS 

In this paper we have investigated whether the sce- 
nario of aging to hydrodynamics depicted in Fig. 1 for 
conventional gases applies to granular gases as well, even 
at high dissipation. Here the term hydrodynamics means 
that the velocity distribution function, and, hence, the 
irreversible fluxes, is a functional of the hydrodynamic 
fields (density, flow velocity, and granular temperature) 
and, thus, it is not limited to the NS regime. To ad- 
dress the problem, we have restricted ourselves to un- 
steady states in two classes of flows, the USF and the 
ULF (see Figs. 2 and 3). While the USF is an incom- 
pressible flow (V • u = 0) and the ULF is a compress- 
ible one (V ■ u ^ 0), they share some physical features 
(uniform density, temperature, and rate of strain ten- 
sor) that allow for a unified theoretical framework. Both 
flows admit heating states (where viscous heating prevails 
over inelastic cooling) and cooling states (where inelastic 
cooling overcomes viscous heating), until a steady state 
is eventually reached. Moreover, in the ULF with posi- 
tive longitudinal rates only "super-cooling" states (where 
inelastic cooling adds to "viscous cooling") arc possible 
and, thus, no steady state exists. 

Two complementary routes have been adopted. First, 



the BGK-likc model kinetic equation (3.1) has been used 
in lieu of the true inelastic Boltzmann equation, which 
allows one to derive a closed set of nonlinear first-order 
differential equations [cf. Eq. (3.8)] for the temporal evo- 
lution of the elements of the pressure tensor. A nu- 
merical solution with appropriate initial conditions and 
elimination of time between the reduced pressure tensor 
P*l and the reduced rate of strain a* provides the non- 
Newtonian functions P*j(a*), from which one can con- 
struct the viscosity function n*(a*) [cf. Eq. (2.28)] and 
(only in the USF case) the viscometric functions ^l(a*) 
[cf. Eq. (2.32)] and #|(a*) [cf. Eq. (2.33)]. The numerical 
task can be avoided at the cost of introducing approxima- 
tions. The one we have worked out consists of expanding 
the solution in powers of a parameter q measuring the 
hardness of the interaction (q = for inelastic Maxwell 
particles and q = \ for inelastic hard spheres), truncat- 
ing the expansion to first order, and then constructing 
Padc approximants. This yields explicit expressions for 
the rate of strain dependence of the rheological functions. 
In the USF case the viscosity and the first viscometric 
functions are given by Eqs. (3.51) and (3.53), respec- 
tively, complemented by Eqs. (3.35), (3.44), and (3.45); 
the second viscometric function vanishes in the BGK-likc 
kinetic model (3.1). As for the ULF, only one rheological 
function (viscosity) exists and it is given by Eq. (3.50), 
complemented by Eqs. (3.36) and (3.48). While one could 
improve the approximation by including terms in the q 
expansion beyond the first-order one [66], the approx- 
imation considered in this paper represents a balanced 
compromise between simplicity and accuracy (see Figs. 4 
and 5). In fact, the results predicted by the kinetic model 
(3.1) and the simplified rheological model described by 
Eqs. (3.50), (3.51), and (3.53) share the non-Newtonian 
solution for q = as well as the steady-state values and 
the values at a* = for arbitrary q. 

The second, and most important, route has been the 
numerical solution of the true Boltzmann equation by the 
stochastic DSMC method for three values of the coeffi- 
cient of restitution and a wide sample of initial condi- 
tions. The most relevant results are summarized in Figs. 
10 (USF) and 14 (ULF). Those figures show an excel- 
lent degree of overlapping of the rheological functions ob- 
tained by starting from the different initial conditions, al- 
though the USF viscometric functions may exhibit large 
fluctuations. The overlapping takes place after a kinetic 
stage lasting about 5 collisions per particle for the heat- 
ing states considered and between 5 and 50 collisions per 
particle for the cooling states considered (see Table III). 
In the latter states the whole dynamics is slowed down 
with respect to the heating states, so the increase of the 
duration of the kinetic stage is correlated with a similar 
increase of the duration of the subsequent hydrodynamic 
stage. We have also observed that the characteristic time 
periods increase as the inelasticity decreases. Figures 10 
and 14 also show that the rheological model, despite its 
simplicity, captures reasonably well, even at a quantita- 
tive level, the main features of the DSMC results. An 
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exception is the USF second viscometric function which, 
although about 10 times smaller than the first viscomet- 
ric function, is unambiguously nonzero. 

In summary, we believe that our results provide strong 
extra support to the validity of a hydrodynamic descrip- 
tion of granular gases outside the quasielastic limit and 
the NS regime. Of course, it is important to bear in 
mind that the USF and ULF are special classes of flows 
where no density or thermal gradients exist, except dur- 
ing the early kinetic stage (see, for instance, Fig. 11). 
Thus, the results presented here do not guarantee a pri- 
ori the applicability of a non-Newtonian hydrodynamic 
approach for a general situation in the presence of den- 
sity and thermal gradients. On the other hand, recent 



investigations for Couette-Fourier flows [18-20, 82, 83] 
nicely complement the study presented in this paper in 
favor of a (non-Newtonian) hydrodynamic treatment of 
granular gases. 
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